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CHAPTER  1 


INTRODUCTION 

The  current  trend  in  microwave  and  millimeter- wave  integrated  circuits  (MICs)  is 
towards  higher  operating  frequencies,  higher  packing  densities,  and  more  stringent 
performance  requirements.  As  a  result,  earlier  models  based  on  approximate  analysis 
techniques  are  becoming  insufficiently  accurate  for  use  in  computer-aided  design  (CAD) 
packages.  In  particular,  simplifying  assumptions,  such  as  the  quasi-static  approximation 
and  the  magnetic  wall  approximation,  are  becoming  inappropriate  in  the  analysis  of 
microstrip  discontinuities  such  as  bends,  steps,  and  stubs.  Instead,  a  rigorous  full-wave 
analysis  is  needed  to  obtain  a  more  accurate  characterization  of  these  passive  circuits.  Such 
an  analysis  captures  the  increasingly  significant  effects  of  coupling,  dispersion,  and 
radiation. 

One  of  the  full-wave  analysis  techniques  that  has  received  considerable  attention 
recently  is  the  spectral  domain  approach.  This  is  an  integral  equation  approach  that  is 
solved  numerically  in  the  spectral  domain  by  applying  the  method  of  moments.  The  open- 
end  and  gap  discontinuities  have  been  analyzed  with  this  technique  in  [1]  and  [2],  and  the 
analysis  of  the  step  junction  has  been  treated  in  [3].  In  a  recent  work  by  Jackson  [4],  the 
analyses  of  the  step,  stub,  and  bent  stub  were  presented. 

Although  the  spectral  domain  approach  has  been  successfully  employed  in  a  variety 
of  problems,  there  is  a  definite  need  for  improvement.  This  approach  is  computationally 
expensive  for  even  the  simplest  geometries  and  quickly  becomes  intractable  for  larger 
problems.  The  primary  source  of  difficulty  is  the  evaluation  of  the  inner  products  arising 
from  the  moment  method.  In  the  spectral  domain  approach,  these  inner  products  are  two- 
dimensional  integrals  with  infinite  limits  of  integration.  Since  these  integrals  are  too 
complicated  to  be  evaluated  analytically,  we  must  resort  to  numerical  techniques. 
Unfortunately,  performing  the  integrations  numerically  is  complicated  by  the  fact  that  the 


integrands  will  have  one  or  more  singularities  and  are  often  highly  oscillatory.  The  most 
troubling  and  time-consuming  aspect  of  these  integrations,  however,  is  the  fact  that  the 
integrands  decay  quite  slowly. 

The  major  portion  of  this  thesis  is  therefore  devoted  to  the  efficient  evaluation  of  the 
inner  products.  In  particular,  an  asymptotic  acceleration  technique  is  applied  to  the 
numerical  integrations  to  enhance  their  rate  of  convergence.  In  this  technique,  the 
asymptotic  form  of  the  integrand  is  derived  and  then  subtracted  from  the  original  integrand. 
The  resulting  integral  will  then  converge  more  rapidly.  In  order  to  recover  the  original 
inner  product ,  the  asymptotic  form  is  integrated  separately  and  added  back.  An  overall 
speedup  in  the  computation  is  realized  only  if  the  asymptotic  form  can  be  integrated 
efficiently.  In  our  work,  the  integral  of  the  asymptotic  form  has  been  derived  analytically 
in  closed  form.  In  this  way,  the  maximal  speedup  from  the  asymptotic  extraction  is 
achieved,  and  more  importantly,  the  technique  is  exact  in  the  sense  that  no  simplifying 
approximations  were  made. 

In  Chapter  2,  the  general  microstrip  discontinuity  problem  is  formulated  in  the 
spectral  domain.  The  moment  method  is  then  applied  in  Chapter  3,  resulting  in  a  system  of 
equations  with  the  inner  products  described  above.  In  Chapter  4,  the  difficulties  involved 
in  the  computation  of  these  inner  products  is  discussed  in  detail.  In  addition,  symmetry 
relations  are  presented  that  allow  a  reduction  in  the  range  of  integration.  The  asymptotic 

acceleration  technique  is  presented  in  Chapter  5.  In  Chapter  6,  comparisons  are  made 

* 

between  an  algorithm  that  utilizes  the  acceleration  technique  and  one  that  does  not.  Finally, 
conclusions  and  directions  for  future  work  are  given  in  Chapter  7. 


CHAPTER  2 


THEORETICAL  DEVELOPMENT 

With  the  increase  in  the  operating  frequencies  of  microwave  integrated  circuits, 
previous  models  of  passive  circuits  based  on  approximate  analysis  techniques  are  becoming 
insufficiently  accurate.  With  this  in  mind,  a  full-wave  analysis  in  the  spectral  domain  is 
developed  in  this  work  to  derive  more  accurate  models  for  a  variety  of  microstrip 
discontinuities.  The  analysis  developed  herein  is  rigorous,  one  that  incorporates  wave 
phenomena  including  radiation,  dispersion,  coupling,  and  surface  wave  propagation.  In 
this  approach,  the  microstrip  discontinuity  is  characterized  by  solving  for  the  reflection  and 
transmission  coefficients.  From  these  coefficients,  a  suitable  circuit  model  can  be  obtained 
for  use  in  circuit  simulation  programs. 

In  this  chapter,  the  conventional  spatial  domain  analysis  is  presented  first.  The 
electric  field  is  represented  in  terms  of  the  currents  on  the  microstrip  through  a  dyadic 
Green's  function.  Since  the  Green's  function  is  unavailable  in  closed  form  in  the  spatial 
domain,  the  formulation  is  continued  in  the  spectral  domain  wherein  the  Green's  function  is 
known.  After  obtaining  the  general  spectral  domain  equations,  these  equations  are  applied 
to  the  microstrip  discontinuity. 

The  microstrip  discontinuity  is  analyzed  as  a  two-port  network  that  is  excited  by  a 
precomputed  incident  current  density.  This  sinusoidally  varying  current  is  determined  by 
solving  for  the  fundamental  mode  on  a  uniform  microstrip  line.  The  forms  of  the  reflected 
and  transmitted  currents  on  the  input  and  output  microstrip  lines  are  also  assumed  to  be  the 
same  as  for  the  current  on  a  uniform  microstrip  line.  Therefore,  the  reflected  and  the 
transmitted  currents  are  precomputed  to  within  a  multiplicative  constant,  namely,  the 
reflection  or  transmission  coefficient.  In  the  process  of  solving  for  the  unknown 
coefficients,  the  currents  within  the  discontinuity  are  also  obtained. 
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2.1  The  Dyadic  Green's  Function 

The  open  microstrip  structure  considered  is  shown  in  Fig.  2.1.  (Figures  and  tables 
appear  at  the  end  of  each  chapter.)  The  conducting  strip  is  located  on  a  substrate  of  relative 
dielectric  permittivity  £r  and  of  thickness  d.  In  conventional  spatial  domain  analysis,  the 
electric  field  is  represented  in  terms  of  the  currents  on  a  conducting  body  through  a  dyadic 
Green's  function.  Assuming  time  harmonic  fields  with  the  time  convention,  the 
electric  field  on  the  substrate  can  be  represented  as 

E(x,z)  =  JJ*G(x  -  x',z  -  z')  J(x',z')  dx'dz'  (2.1) 


where  J  is  the  current  density,  and  G  is  the  dyadic  Green’s  function.  The  dyadic  Green's 
function  G(x,z)  represents  the  electric  field  at  the  point  (x,z)  due  to  an  infinitesimal  dipole 
of  unit  strength  located  at  the  origin.  The  integration  in  (2.1)  is  performed  over  the  strip 
where  the  current  exists. 

When  Equation  (2.1)  is  enforced  on  the  strip,  the  tangential  electric  field 
components  Ex  and  Ez  must  be  zero  to  satisfy  the  boundary  condition.  This  leads  to  a  set 
of  coupled  homogeneous  integral  equations 


JJ[Gxx(x  -  x'.z  -  z')Jx(x',z') 

+  Gxz(x  -  x',z-  z')Jz(x',z')]  dx'dz'  =  0 


JJ[Gzx(x  -  x',z  -  z')Jx(x',z') 

+  Gzz(x  -  x',z  -  z')Jz(x',z')]  dx'dz'  =  0 


(2.2) 


where  the  vector  equation  has  been  broken  into  its  scalar  components.  If  the  Green’s 
function  were  known  in  the  spatial  domain.  Equation  (2.2)  could  be  used  to  solve  for  the 


currents  on  the  microstrip.  Unfortunately,  for  the  microstrip  structure,  a  closed-form 
expression  for  the  Green's  function  is  not  known. 

A  parallel  analysis  may  be  carried  out  in  the  spectral  domain  where  the  Green’s 
function  for  a  grounded  dielectric  slab  can  be  obtained  using  a  transmission  line  approach 
[5].  Taking  the  Fourier  transform  of  (2.1)  converts  the  convolution-type  equations  to  the 
following  algebraic  equations 

Gxx(a,(3)jx(a,p)  +  Gxz(a,p)Iz(a,|3)  =  Ex(a,p) 

G2X(a,P)Jx(a,13)  +  Gzz(a,p)Jz(a,p)  =  Ez«x,(3)  v 


where  the  tilde  denotes  the  Fourier  transform  operator.  The  definition  of  the  Fourier 
transform  and  its  inverse  transform  have  been  taken  to  be 


and 


f(oc,p)=  J  J  f(x,z)e~j(ax+3z)  dxdz  (2.4) 


f(X,Z) 


—Ur  f  f  f(a,P)ej(ca+pz)  dadp, 

(2k)  _UjL 


(2.5) 


respectively. 

Note  that  Equation  (2.1)  was  transformed  instead  of  Equation  (2.2).  This  is 
necessary  since  the  transform  is  taken  over  the  entire  spatial  domain  and  (2.2)  is  defined 
only  on  the  strip.  As  a  result,  the  transform  of  the  electric  field  in  (2.3)  is  nonzero.  This 
electric  field  is  unknown,  but  can  be  eliminated  when  the  moment  method  is  applied  (this  is 
discussed  further  in  Chapter  3). 

The  spectral  domain  Green's  function  for  the  open  microstrip  structure  is  derived  in 
detail  by  Itoh  and  Menzel  [6];  therefore,  only  the  result  is  presented  here.  In  the  derivation, 
the  strip  is  considered  to  be  infinitesimally  thin,  and  both  the  ground  plane  and  the  strip  are 


treated  as  perfect  electric  conductors.  An  equivalent,  more  compact  form  of  the  Green's 
function  is  given  by  [7] 
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GxzCa.P) 


a2D3  -k2D} 

£oeoDiD2 

apP3 

G)e0D1D2 


Gzx(a,p)  =  Gxz(a,p) 


Gzz(a,p) 


P2D3-feoD| 

CO£0D1D2 


(2.6) 


(2.7) 

(2.8) 
(2.9) 


where 


s  =  3/k2-a2-P2 
s'  =  ^er k2  -  a2  -p2 

Dj  =  s'- jsercot(s'd) 
P2  =  s- js'cot(s'd) 
D3  =  s' -  jscot(s'd) 


and  where  the  positive  real  or  negative  imaginary  root  is  taken  in  the  evaluation  of  s  and  s'. 

2.2  Spectral  Equations  Describing  the  Microstrip  Discontinuity 

The  geometry  of  a  general  microstrip  discontinuity  is  shown  in  Fig.  2.2,  where  the 
particular  discontinuity  is  defined  within  a  discrete  mesh  region.  The  discontinuity  is 
analyzed  as  a  two-port  network  where  the  reference  planes  are  defined  at  z=0  and  z=L. 
The  network  is  excited  at  z=0  by  an  incident  current  density  which  is  assumed  to  have  the 
same  form  as  the  current  density  on  a  uniform  microstrip  line.  The  reference  planes  are 


assumed  to  be  sufficiently  far  from  the  discontinuity  so  that  the  reflected  and  the  transmitted 
currents  can  be  approximated  by  currents  on  a  uniform  microstrip  line.  The  currents  within 
the  mesh  region  are  entirely  unknown.  The  spectral  current  density  appearing  in  Equation 
(2.3)  is  then  the  superposition  of  the  boundary  currents  and  the  mesh  region  currents, 
given  by 


J(a,p)  =  JW)  +  rjR«x,p)  +  TjT(a,P)  +  jmr(a,p)  (2. 10) 

T  R  X 

where  J  is  the  incident  current  density,  J  and  J  are  the  reflected  and  transmitted 
currents,  and  Jmr  is  the  unknown  current  density  within  the  mesh  region.  The  reflection 
and  transmission  coefficients  are  given  by  T  and  T,  respectively. 

The  boundary  current  densities  in  (2.10)  may  be  represented  in  the  spatial  domain 
as 


J1  (x, z)  =  [x  fxi(x,z)  +  z  fzl(x,z)]  e-JP°lZ 

(2.11) 

JR(x,z)  =  [-x  fxl(x,z)  +  z  fzl(x,z)]  e+JP°lZ 

(2.12) 

JT(x,z)  =[x  fx2(x,z)  +  z  fz2(x,z)]  e~jPo2Z 

(2.13) 

where  the  functions  fxi,fzl.fx2.  and  fz2  represent  the  transverse  distribution  of  currents, 
and  (3oi  and  po2  are  the  propagation  constants.  The  propagation  constants  and  the 
transverse  distributions  are  precomputed  by  solving  for  the  fundamental  mode  of  a  uniform 
microstrip  line  [5].  Specifically,  Equation  (2.3)  is  solved  numerically  for  the  case  in  which 
the  microstrip  is  invariant  in  the  z-direction  and  there  is  no  source  excitation.  The 
numerical  technique  is  essentially  a  simplification  of  the  technique  to  be  described  in 
Chapter  3. 

The  only  remaining  unknowns  in  (2.10)  are  Jmr,  T  and  T.  These  are  found  by 
solving  (2.3).  Substitution  of  (2.10)  into  (2.3)  yields 
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GxJr+G-2jf->-r(G5>J*  +G«j? 

dJT +GJ™  +  r(a„3?  +o22i? 


j  +  T|GXXJX  +  Gx2J2 
)  +  t(G2xJ)+G2zJJ 


)  "  '  Grjl,  ~G,>; I,  j 

)  =  e2-(g2XJ‘x+g22J‘) 


(2.14) 


where  for  simplicity  the  arguments  of  the  functions  have  been  dropped.  The  numerical 
solution  of  (2.14)  is  discussed  in  Chapter  3. 


y 


ground  plane 


Figure  2.1.  Cross  section  of  an  open  microstrip  line. 


Figure  2.2.  Geometry  of  a  general  microstrip  discontinuity. 
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CHAPTER  3 

NUMERICAL  PROCEDURE 

In  Chapter  2,  the  general  microstrip  discontinuity  problem  was  formulated  in  the 
spectral  domain,  and  the  electric  field  was  represented  in  terms  of  the  current  densities  on 
the  strip.  In  this  chapter,  the  moment  method  is  applied  to  the  spectral  domain  Equation 
(2.14)  to  generate  a  system  of  equations  that  can  be  solved  numerically.  The  first  section 
of  this  chapter  introduces  the  moment  method  procedure.  In  the  following  section,  the 
procedure  is  applied  to  the  discontinuity  problem,  and  the  current  densities  in  (2.14)  are 
expanded  in  terms  of  basis  functions.  The  final  section  provides  the  criteria  for  selecting  a 
suitable  set  of  basis  functions  and  presents  the  basis  functions  used  in  this  work. 

3.1  Overview  of  the  Moment  Method 

The  moment  method  is  a  mathematical  technique  used  to  reduce  functional 
equations  to  matrix  equations  [8].  Consider  the  linear  equation 

L(f)*s  (3.1) 

In  this  equation,  L  is  a  linear  operator,  s  represents  a  source,  and  f  is  an  unknown  function 
such  as  a  field  or  current.  The  problem  is  to  determine  the  function  f  when  the  inverse 
linear  operator  L'1  is  unknown. 

The  function  f  is  first  expanded  into  a  series  of  basis  functions  fj.fo . f\r  defined 

over  the  domain  of  L 


N 

f  =  Xaifi 

i=l 


(3.2) 
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The  basis  functions  should  be  linearly  independent  and  chosen  so  that  a  weighted  sum  of 
these  functions  can  represent  f  reasonably  well.  Depending  on  the  choice  of  basis 
functions,  the  series  representation  will  be  approximate  or  exact .  The  basis  functions  used 
in  this  work  are  subsectional  basis  functions  where  each  fj  is  defined  over  only  a 
subsection  of  the  domain  of  f.  In  this  case,  Equation  (3.2)  will  most  likely  be  an 
approximate  representation  of  f,  given  a  finite  N. 

Equation  (3.2)  is  then  substituted  into  the  Linear  Equation  (3.1)  to  give 

N 

2>iL(fi)  =  s  (3.3) 

i=l 


A  set  of  N  equations  is  needed  to  solve  for  the  unknown  coefficients  These 

equations  are  obtained  by  choosing  N  weighting  or  testing  functions  w^wt . wN  and 

taking  a  suitable  inner  product  of  (3.3)  with  each  testing  function.  When  the  basis  and 
testing  functions  are  the  same,  the  procedure  is  known  as  Galerkin's  method.  The  final  set 
of  equations  is  expressed  in  matrix  form  as 

(w^Lffj))  (w1,L(f2))  •••  (wltL(fN)) 

<w2,L(f1))  (w2,L(f2))  •••  (w2,L(fN)} 

_(wN’L(fi))  (w^,L(f2)}  •••  (wjvf.Lff^)) 

These  equations  can  be  solved  numerically  using  a  number  of  techniques.  Since  the 
limiting  factor  in  this  method  is  the  computation  time  required  to  evaluate  the  inner 
products,  the  size  of  the  matrix  in  (3.4)  must  be  kept  small.  Consequently,  the  matrix 
inversion  is  not  a  problem. 


*  al  * 

(W1’S)' 

a2 

= 

(w2,s) 

_aN_ 

_{wN,s) 

(3.4) 


12 


3.2  Application  of  the  Moment  Method  to  the  Microstrip  Discontinuity 

The  moment  method  is  now  applied  to  the  spectral  domain  Equation  (2.14) 
describing  the  microstrip  problem.  Here,  the  Green's  function  is  the  linear  operator  and  the 
unknown  function  is  the  current  density  in  the  discontinuity  region.  Two  additional 
unknowns  are  the  reflection  and  transmission  coefficients  T  and  T.  The  source  function  is 
the  incident  current  density  J1. 

To  apply  the  moment  method,  the  spectral  current  densities  in  the  mesh  region  are 
first  expanded  in  terms  of  basis  functions 

N 

J™r(a,B)  =  XaiJxi(<xP) 

M  <3-5> 

If(a,p)  =  XbiJzi(a,P) 

i=l 


where  a;  and  bi  are  the  weight  coefficients  and  Jxj  and  J2j  are  the  basis  functions.  The 
boundary  currents  J1,  JR,  and  JT  are  also  expanded  in  terms  of  basis  functions  such  that 


]J(a,P)  =  Xc?j?i(a,P) 

1  (3.6) 

I?(a,P)  =  ]Td?ig(a,P) 
i 

where 

Q  =  I,R,  or  T 

All  of  the  constants  c^  and  d^  are  known,  since  the  boundary  currents  are  precomputed. 
The  boundary  current  densities  are  actually  represented  in  terms  of  basis  functions  for 
reasons  that  are  not  associated  with  the  moment  method.  The  main  reason  is  that  the 
transverse  current  distribution  is  calculated  numerically  and  therefore  must  be 


approximated.  Also,  the  boundary  currents  are  semiinfinite,  but  they  are  truncated  to 
simplify  the  computation. 

The  current  expansions  (3.5)  and  (3.6)  are  then  substituted  into  (2.14),  and  the 
inner  product  of  this  equation  is  taken  with  testing  functions.  The  inner  product  definition 
adopted  in  this  work  is 


(f(a,P),g(cc,p))  =  J  J  f*(a,P)g(a,p)  dadfi  (3.7) 


where  the  function  f  is  the  testing  function,  and  the  function  g  is  the  Green's  function 
multiplied  by  a  current  basis  function.  Physically  this  type  of  testing  will  give  the 
integrated  electric  field  over  the  area  of  the  testing  function  due  to  a  current  element 
radiating  over  the  grounded  dielectric  slab. 

The  testing  procedure  is  almost  Galerkin.  The  first  row  in  (2.14)  is  tested  N  times 
with  Jxi(a,p),  and  the  second  row  is  tested  M  times  with  Jzj(a,P).  However,  this  will 
yield  only  (N+M)  equations  to  solve  for  the  (N+M+2)  unknowns.  The  last  two  equations 
are  obtained  by  testing  at  the  junctions  between  the  mesh  region  and  the  input  and  output 
reference  planes.  Since  the  z-directed  current  will  be  dominant,  the  inner  product  testing  is 
taken  with  the  second  row  of  (2.14).  The  boundary  testing  functions  are  denoted  as 
Wj(a,p)  and  W2(a,p).  The  final  matrix  form  is  then 


ZXX  ZXZ 

ZxR  ZxT 

V 

r  vxii 

ZZX  ZZZ 

ZzR  ZzT 

bi 

VZI 

Zwlx  Zwlz 

ZwlR  ZwlT 

r 

VwlI 

-Zw2x  Zw2z 

Zw2R  Zw2T- 

_T. 

-Vw2I- 

where  the  matrix  elements  shown  are  submatrices. 
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The  upper  left  block  of  submatrices  in  (3.8)  results  from  the  testing  between  basis 
functions  within  the  mesh  region.  These  are  the  largest  submatrices  and  increase  in  size  as 
the  square  of  the  number  of  basis  functions.  The  upper  right  block  consists  of  inner 
products  where  the  testing  is  between  the  mesh  region  basis  functions  and  the  reflected  or 
transmitted  currents.  Each  of  these  submatrices  are  column  vectors.  The  last  two  rows  in 
(3.8)  are  the  two  additional  equations  needed  to  solve  for  T  and  T.  Finally,  the  vector  on 
the  right-hand  side  results  from  taking  the  inner  product  between  each  of  the  testing 
functions  and  the  incident  current. 

The  similarities  between  different  submatrices  are  readily  seen  when  the 
submatrices  are  written  in  compact  form.  The  general  term  in  each  submatrix  is  given  by 


where 


(Zpq)ji  =  (jpj’GpqJqi) 

(matrix) 

(ZpQ)j  =  (jpj»(GpXJx  +<-*pzj?  )) 

(column  vector) 

(Zwnq)i  =:  (^Yn’^zqjqi) 

(row  vector) 

(Zwi,q)  =  ^Wn,^GpxJ^  +GpZj^^ 

(scalar) 

(Vpl)j=(jpj’(Gpxjx+Gpzj5)) 

(column  vector) 

(^wnl)  =  /wn,^GzxJx  +  GZZJZ  j\ 

(scalar) 

p  =  x  or  z 
q  =  x  or  z 
Q  =  RorT 
n  =  1  or  2 


Note  that  the  unknown  electric  field  in  the  right-hand  side  of  (2.14)  was  eliminated 
in  the  testing  process  and  does  not  appear  in  the  blocks  Vpi  or  VWnI-  This  is  explained  by 
applying  Parseval's  theorem  to  the  inner  products.  As  an  example,  consider  the  inner 


product  <Jxj,GxxJxi>.  which  is  equivalent  to  <Jxj,Ex>.  After  applying  Parseval's 
theorem,  the  inner  product  becomes 
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(ixj.GxxJxi)  =  J  pxj(a,p)Ex(cc,p)dadp 

— oo  — oo 

oo  oo 

=  4j i2  J  Jlxj(x,z)Ex(x,z)  dxdz  (3.9) 

—  oo  —  oo 

=  0 

The  inner  product  is  zero,  because  in  the  spatial  domain  the  current  density  and  the  electric 
field  on  the  substrate  are  zero  in  complementary  regions.  For  instance,  the  current  is  zero 
everywhere  except  on  the  conducting  strip  where  the  electric  field  is  zero.  Therefore,  the 
integrand  in  the  second  integral  is  identically  zero. 

3.3  Basis  Functions 

As  mentioned  in  Section  3.1,  the  choice  of  basis  functions  can  affect  the  accuracy 
of  the  current  representation.  A  good  choice  can  also  reduce  the  number  of  basis  functions 
needed  to  obtain  adequate  results.  When  choosing  basis  functions,  there  are  several  issues 
that  must  be  considered.  These  are  the  geometry  of  the  problem,  the  current  behavior,  and 
the  inner  product  evaluation. 

The  geometry  of  the  problem  will  dictate  the  use  of  either  entire  domain  or 
subdomain  basis  functions.  In  this  work,  subdomain  basis  functions  are  used  so  that  a 
variety  of  microstrip  discontinuities  may  be  analyzed.  The  geometry  of  the  discontinuity  is 
actually  defined  by  the  location  of  the  current  basis  functions,  since  the  basis  functions  are 
placed  only  on  the  conducting  strip  where  current  exists.  For  complicated  geometries,  the 
subdomain  basis  functions  also  offer  the  advantage  that  no  prior  knowledge  of  the  current 
behavior  is  needed.  This  is  contrasted  with  the  patch  antenna  problem  where  entire  domain 


basis  functions  may  be  used,  because  the  geometry  is  fixed  and  there  is  some  knowledge  of 
the  modal  distribution  of  current. 
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In  order  to  model  the  currents  efficiently,  some  important  aspects  of  the  current 
behavior  should  be  built  into  the  basis  functions.  One  such  characteristic  is  the  well-known 
edge  condition.  This  refers  to  the  fact  that  for  current  traveling  along  an  edge,  the  current 
distribution  will  tend  towards  infinity  near  the  edge.  Another  characteristic  is  that  currents 
traveling  towards  an  edge  of  a  conductor  should  go  to  zero  by  current  continuity.  The 
basis  functions  are  chosen  to  exhibit  these  properties. 

The  final  consideration  in  choosing  the  basis  functions  is  to  ensure  that  the 
integration  in  the  inner  product  evaluation  is  convergent.  A  careful  look  at  Equations  (2.6)- 
(2.9)  will  show  that  the  Green's  function  increases  without  bound  for  large  a  and  p.  If  the 
two-dimensional  integration  is  convened  into  polar  coordinates  by  making  the  substitutions 
a=psin0  and  p=pcos0,  it  can  be  shown  that  for  large  p  the  components  of  the  asymptotic 
Green’s  function  tend  to  Gxx  =  psin20,  Gxz  =»psin0cos0,  and  Gzz  *  pcos20.  Since  the 
coordinate  transformation  introduces  another  p  as  a  metric  coefficient,  the  spectral  domain 
basis  and  testing  functions  must  introduce  a  factor  of  at  least  p'p,  p>3,  to  guarantee 
convergence. 

In  the  spatial  domain,  the  "rooftop"  expansion  functions  depicted  in  Fig.  3.1  are 
used  as  basis  functions  in  the  mesh  region.  These  basis  functions  are  composed  of  a 
triangle  function  in  the  direction  of  current  flow  and  a  pulse  function  in  the  other  direction. 
The  pulse  and  triangle  functions  are  defined  as 


nT(0  = 


(3.10) 
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^w(0  — 


i-M, 


w 


0, 


|t|<w 
|t|>  w 


(3.11) 


where  the  corresponding  Fourier  transforms  with  respect  to  the  spectral  variable  co  are 


n-r(co)  =  T--n-^-T— ) 
1  (coT/2) 


Aw(co)  =  w 


sin  (cow  12) 
(cow  /  2) 


i2 


(3.12) 

(3.13) 


Combining  these  two  functions  in  the  spatial  domain,  the  mesh  region  basis 
functions  for  the  x  and  z-directed  currents  become 


Jxi(x,z)  =  AWx(x-xi)ntx(z-zi)  (3.14) 

Jzi(x,z)  =  ntz(x-xi)AWz(z-zi)  (3.15) 


For  each  component  of  current,  the  widths  and  thicknesses  of  all  the  "rooftop"  basis 
functions  are  assumed  to  be  equal.  The  center  of  the  ith  basis  function  is  (xi,z,),  and  the 
basis  functions  are  positioned  so  that  the  triangle  functions  of  neighboring  "rooftops" 
overlap  and  the  pulse  functions  are  edge  to  edge.  The  current  represented  by  these  basis 
functions  goes  to  zero  in  the  direction  of  current  flow,  and  the  basis  functions  are  suited  to 
model  the  edge  condition. 

The  spectral  domain  basis  functions  are  obtained  by  taking  the  Fourier  transform  of 
(3.14)  and  (3.15) 


Jxi(a,p)  =  AWx(a)fltx(P)  e-j(aXi+Pzi)  (3.16) 

Jzi(a,P)  =  fltz(a)AWz(p)e-j(ax>+Pzi) 


(3.17) 
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Each  of  these  basis  functions  introduces  a  p'3  term  in  polar  coordinates;  therefore,  the  inner 
product  integrals  are  convergent  when  the  testing  procedure  is  Galerkin. 

The  boundary  currents  are  modeled  in  a  manner  consistent  with  the  mesh  region 
basis  functions.  The  transverse  current  distribution  is  approximated  as  a  weighted  sum  of 
pulse  functions  for  the  z-directed  currents  and  triangle  functions  for  the  x-directed  currents. 
The  weighting  coefficients  are  found  by  sampling  the  precomputed  current  distribution  at 
the  center  of  each  basis  function.  The  x-directed  current  is  truncated  at  the  mesh  region  as  a 
pulse  function  while  the  z-directed  current  terminates  in  the  mesh  region  as  a  triangle 
function.  These  basis  functions  are  multiplied  by  the  propagation  factor  e+iPoZ  and  are 
truncated  at  a  distance  tb  from  the  mesh.  The  boundary  currents  in  the  spatial  domain  are 
represented  as 


where 


J?(x,z)  =  Xc?  Awx(x  "  x0ntb(z  +  /  2)  e+Jp0lZ 

JQ(x,  z)  =  ntz(x  -  Xj)P(z  +  tb  /  2)  e+jPolZ 
i 

jJ (x,z)  =  Awx(x  -  xi)ntb(z -  L  -  tb  /  2)  e-jPo2Z 
i 

J I (x,z)  =  £d?  ntz(x  -  Xj)p(z -  L  -  tb  /  2)  e_JPo2Z 


Q  =  I  or  R 


and  where 


P(t)  = 


— [(tb/2  +  wz)-|t|  ], 


w 


1, 


tb/2  <|t|  c(tb/2  +  wz) 
|t|^tb/2 


(3.18) 


(3.19) 


In  (3.18)  the  negative  sign  is  taken  for  the  incident  current  and  the  positive  sign  for  the 
reflected  current.  Figures  3. 2-3.5  illustrate  how  the  input  boundary  currents  J1  and  JR  are 
represented  in  terms  of  these  basis  functions. 
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The  boundary  testing  functions  are  also  chosen  to  be  compatible  with  both  the  mesh 
region  basis  functions  and  the  boundary  current  basis  functions.  These  testing  functions 
are  a  weighted  sum  of  z-directed  "rooftop"  basis  functions,  where  the  weights  are  chosen 
to  be  consistent  with  the  z-directed  boundary  currents.  The  testing  function  Wi  straddles 
the  z=0  reference  plane  at  the  junction  between  the  mesh  region  and  the  boundary  currents. 
Similarly,  the  testing  function  W2  straddles  the  reference  plane  at  z=L.  These  boundary 
testing  functions  are  represented  in  the  spatial  domain  as 

Wi(x,z)  =  £d[  ntz(x-Xi)AWz(z)  (3.20) 

i 

W2(x,z)  =  nt2(x-Xi)AWz(z-L)  (3.21) 

i 

Finally,  the  basis  and  testing  functions  presented  in  this  section  are  combined  with 
the  moment  method  procedure  in  Section  3.2  to  construct  the  set  of  Equations  (3.8).  The 
numerical  evaluation  of  the  inner  products  is  discussed  next  in  Chapter  4. 
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Figure  3.2.  Transverse  distribution  of  z-directed  boundary  currents  represented  by  8  basis 
functions. 


✓ 


mesh  region 

I  R 

Figure  3.3.  Representation  of  Jz  or  Jz  with  3  basis  functions. 


Figure  3.4.  Transverse  distribution  of  x-directed  boundary  currents  represented  by  7  basis 
functions. 


mesh  region 


I  R 

Figure  3.5.  Representation  of  Jx  or  Jx  with  3  basis  functions. 
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CHAPTER  4 

INNER  PRODUCT  EVALUATION 

The  formulation  of  the  microstrip  discontinuity  problem  is  quite  elegant  and  simple; 
implementing  the  numerical  procedure  is  considerably  more  difficult.  At  first  glance, 
constructing  the  set  of  matrix  equations  in  (3.8)  appears  relatively  straightforward.  The 
only  requirement  is  to  compute  the  inner  products  between  a  testing  function  and  the 
Green's  function  multiplied  by  a  basis  function.  However,  it  is  the  evaluation  of  these 
inner  products  that  is  the  primary  source  of  difficulty  in  applying  this  method. 

There  are  certain  problems  associated  with  the  inner  product  evaluation  that  are 
apparent  immediately.  For  instance,  the  required  integrals  are  too  complicated  to  be 
evaluated  analytically;  therefore,  the  integrations  must  be  performed  numerically.  These 
integrals  are  doubly  infinite  which  creates  an  additional  problem.  Since  the  infinite 
integrations  must  eventually  be  truncated,  the  rate  of  convergence  is  extremely  important. 
Another  problem  is  that  the  Green's  function  introduces  singularities  in  the  integrand, 
further  comp  lie- ting  the  integration. 

With  these  problems  in  mind,  an  algorithm  was  developed  to  perform  the  analysis 
described  in  Chapters  2  and  3.  This  program  was  written  to  construct  the  matrix  equations 
in  (3.8)  efficiently.  All  of  the  inner  products  were  computed  in  parallel,  reducing  the 
number  of  times  the  functions  were  evaluated.  Essentially,  the  program  integrated  the 
entire  matrix  and  column  vector  in  (3.8)  at  the  same  time.  In  the  algorithm,  the 
convergence  problem  was  handled  by  choosing  higher-order  basis  functions  than  those 
described  in  Chapter  3  in  order  to  increase  the  rate  of  decay  in  the  spectral  domain.  After 
running  the  program,  it  was  found  that  the  rate  of  convergence  was  insufficient.  It  wa: 
also  discovered  that  the  singularity  in  the  integrand  was  treated  unsuccessfully  by  adding 
loss  to  the  permittivity  of  the  dielectric  substrate. 


These  findings  warrant  a  more  careful  look  into  the  evaluation  of  the  inner 
products.  Since  the  sizes  of  the  four  upper  left  blocks  in  (3.8)  increase  as  the  number  of 
basis  functions  squared,  these  inner  products  are  considered  in  particular.  In  this  chapter, 
the  symmetry  in  the  integrand  is  examined  and  exploited  to  reduce  the  overall  computation 
time.  Then  the  nature  of  the  singularities  and  possible  methods  to  integrate  them  accurately 
are  discussed.  Finally,  the  convergence  rate  of  the  inner  product  integrals  is  discussed. 

4.1  Symmetry  in  the  Integrands 

Since  the  number  of  inner  products  that  need  to  be  evaluated  in  a  particular  problem 
may  become  large,  it  is  important  to  study  methods  to  compute  the  inner  product  integrals 
more  efficiendy.  If  the  symmetry  in  the  integrands  is  taken  into  account,  the  numerical 
computation  time  of  each  inner  product  can  be  reduced  by  a  factor  of  4.  This  is 
accomplished  by  changing  the  limits  of  integration  from  the  entire  (a,(3)  plane  to  only  the 
first  quadrant.  In  this  work,  the  numerical  integrations  are  actually  performed  in  the  polar 
coordinate  system,  since  the  doubly  infinite  limits  are  then  reduced  to  an  infinite  limit  on 
only  the  p  variable. 

Although  the  integrations  are  performed  in  polar  coordinates,  it  is  easier  to 
investigate  the  symmetry  of  the  integrand  in  the  Cartesian  coordinate  system.  The 
symmetry  of  each  function  in  the  integrand  is  first  analyzed,  and  then  the  entire  integrand  is 
considered.  The  Green's  function  (2.6)-(2.9)  exhibits  the  following  symmetry: 

Gxx(a,(3)  =  Gxx(-a,p)  =  Gxx(-a,-p)  =  G^c^-p) 

Gxz(a,p)  =  -Gxz(-a,p)  =  Gxz(-a,-p)  =  -Gxz(a,-p)  (4.1) 

Gzz(a,p)  =  GqHx.P)  =  GK(-a,-P)  =  Gzz(a,-p) 


The  components  Gxx  and  Gzz  are  symmetric  in  all  four  quadrants,  while  Gxz  is  asymmetric 
in  the  2nd  and  4th  quadrants. 
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The  symmetry  of  the  basis  and  testing  functions  is  determined  by  examining  each 
function  in  Equations  (3.16)  and  (3.17).  The  triangle  and  pulse  functions  are  even 
functions  with  respect  to  a  and  p.  The  function  that  complicates  the  symmetry  is  the 
complex  exponential  that  results  from  the  shift  of  the  basis  functions  from  the  origin  in  the 
spatial  domain.  For  simplicity,  the  exponentials  from  the  testing  function  and  the  basis 
function  are  grouped  together  and  denoted  by 

4 

Z(a,P)  =  e--i(ax+Pz)  (4.2) 

where 

x'  -  (x,  -Xj) 
z'  =  (Zi  -Zj) 

The  location  of  the  testing  function,  or  field  point,  is  (xj.zp  and  the  location  of  the  current 
basis  function,  or  source  point,  is  (x^Zj).  Note  that  the  change  of  sign  on  the  Xj  and  Zj  of 
the  testing  function  results  from  the  conjugation  in  the  inner  product  definition.  The 
complex  exponential  function  has  the  following  symmetry 


Z(-a,-P)  =  Z*(a,P)  (43) 

Z(a,-p)  =  Z*(-a,P) 

Using  the  symmetries  in  (4.1)  and  (4.3),  the  inner  product  <Jxj>GxxJxi>  can  be 
written  as 


(Jxj,GxxJxi)  =  |  ||A^x(a)n?x(p)Gxx(a,P)[z(a,p)  +  Z*(a,p) 

00 


+  Z(-a,P)  +  Z*(-a,P)  jdadp 


(4.4) 
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where  the  integration  is  now  over  only  the  first  quadrant.  Each  of  the  exponential  terms 
corresponds  to  one  quadrant  in  the  (a,P)  plane.  These  terms  can  be  combined  by  grouping 
complex  conjugate  pairs  to  obtain  twice  the  real  part.  The  resulting  cosine  functions  are 
then  combined  using  a  trigonometric  identity.  The  final  form  of  the  inner  product  is  then 

oooo 

(JXj,Gxxjxi)  =  4  J"  J Awx (a) n?x (P) Gxx (a, (3) cos(a x') cos(P z')  dadp  (4.5) 
0  0 


The  limits  of  integration  for  the  other  inner  products  are  reduced  to  one  quadrant  in 
a  similar  fashion  and  are  given  by 


(jxj,GxzIzi)  =  -4  JjAWx(a)fTtx(P)ntz(a)AW2(p) 

00 

x  G^Ca.PjsinfaxOsinCPz')  dadp  (4.6) 

oooo 

(jzj.GzzJzi)  =  4  J  JnJz(a)A2Wz(p)Gzz(a,p)cos(ax,)cos(pz')  dadp  (4.7) 
0  0 


Since  the  components  and  Gzx  of  the  Green's  function  are  identical,  the  inner  product 
<J2j,GzxJxi>  is  also  given  by  Equation  (4.7). 

As  previously  mentioned,  all  of  the  numerical  integrations  arc  performed  in  polar 
coordinates.  In  this  case,  the  limits  of  the  integration  are  (0,7t/2)  in  0  and  (0,°°)  in  p.  One 
final  observation  is  that  when  x'  or  z  is  zero  in  Equation  (4.7),  the  inner  product  is 
identically  zero.  This  agrees  with  the  physical  intuition  gained  from  studying  a  dipole 
radiating  in  free  space.  For  instance,  a  z-directed  dipole  located  at  the  origin  will  produce  a 
radiation  pattern  that  has  no  x-directed  electric  field  along  either  the  x-axis  or  the  z-axis. 
This  property  apparently  remains  true  for  a  dipole  radiating  over  a  grounded  dielectric  slab. 
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4.2  Singularities 

One  of  the  major  difficulties  in  the  evaluation  of  the  inner  products  is  the  numerical 
integration  of  the  singularities.  The  singularities  in  the  integrand  are  introduced  by  the 
terms  Dj  and  D2  in  the  denominator  of  the  Green's  function  (2.6)-(2.9).  Each  of  the  poles 
of  the  Green's  function  corresponds  to  a  resonance  condition,  namely,  a  source-free 
solution  to  Maxwell's  equations  for  the  grounded  dielectric  slab.  Physically,  the  zeros  of 
Dl  and  D2  represent  the  TM  and  TE  surface  wave  poles,  respectively. 

These  surface  wave  poles  are  determined  by  solving  for  the  roots  of  the  equations 

s-  js'cot(s'd)  =  0  (TE)  .  (4.8) 

s'- jsercot(s'd)  =  0  (TM)  (4.9) 

where 

s  =  Vk5"P2 
s' =  Verko~P2 


These  equations  are  functions  of  only  the  variable  p  in  polar  coordinates,  and  therefore  the 
singularities  are  encountered  only  in  the  p  integration. 

A  careful  examination  of  Equations  (4.8)  and  (4.9)  will  show  that  in  the  case  of  a 
lossless  dielectric  (er  real),  the  singularities  occur  for  real  values  of  p  in  the  region 
k0<p<Vei-k0.  The  number  of  singularities  within  this  region  can  be  shown  to  be  [9] 


and 


Nte  = 


[0, 

In, 


X  <  7X/2 

(n  -  l/2)7t  <  x  <  (n  +  l/2)7t,  n  =  1,2,3,- 


NTM  =  {n  +  l,  nrt  <  x  <  (n  +  l)jt,  n=  0,1,2,— 


(4.10) 


(4.11) 


where 


x  =  M^/Er-l 
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Since  the  first  TM  surface  wave  has  a  zero  cutoff  frequency,  there  is  at  least  one  singularity 
in  the  integration.  Additional  surface  wave  poles  are  introduced  as  the  electrical  thickness 
of  the  dielectric  increases.  In  most  applications,  these  higher-order  surface  waves  are 
avoided. 

Since  there  will  always  be  at  least  one  singularity  present,  it  is  necessary  to  develop 
a  suitable  technique  to  numerically  evaluate  the  singular  integrals.  There  are  a  couple  of 
simple  techniques  that  could  be  attempted.  One  is  to  add  loss  to  the  dielectric,  thereby 
shifting  the  poles  slighdy  off  the  path  of  integration.  Another  possibility  is  to  deform  the 
path  of  integration  around  the  singularity  by  applying  the  Cauchy-Goursat  theorem. 
However,  it  will  be  shown  that  the  first  method  is  inadequate,  and  the  second  method  is  not 
theoretically  justified  for  the  inner  product  definition  (3.7). 

The  simplest  technique  to  implement  is  adding  loss  to  the  dielectric  substrate.  The 
dielectric  loss  is  easily  included  in  the  analysis  by  replacing  Ej.  in  (2.6)-(2.9)  by  Ejfl-jtanS). 
where  tan5  is  the  loss  tangent  of  the  dielectric  substrate.  This  loss  will  move  the  poles  off 
the  real  p  axis  so  that  the  poles  will  occur  at  p=p’-jp”,  p">0.  Since  the  path  of  integration 
is  on  the  real  p  axis,  the  integral  is  theoretically  nonsingular. 

For  realistic  problems,  however,  the  integrand  is  still  poorly  behaved  for  values  of 
p  near  p=p’,  because  the  poles  do  not  move  far  off  the  real  p  axis  unless  the  dielectric 
substrate  is  unrealistically  lossy  and/or  thick.  For  example,  consider  a  typical  microstrip 
substrate  of  alumina  (^=9.6,  tan5=10'4)  that  is  25  mils  thick.  An  actual  microstrip  circuit 
designed  on  this  substrate  might  operate  at  frequencies  between  0  and  20  GHz.  For  this 
case,  the  location  of  the  first  TM  surface  wave  pole  is  plotted  in  Fig.  4.1  for  various 
substrate  thicknesses  and  dielectric  losses.  The  substrate  thicknesses  dAo=0.02,  0.04,  and 
0.08  correspond  to  the  operating  frequencies  of  9.5,  18.9,  and  37.8  GHz,  respectively. 
For  each  of  these  thicknesses,  the  real  and  imaginary  pans  of  the  pole  are  shown  for  the 
dielectric  loss  tangents  of  10'6,  10”4  and  10'2.  The  figure  shows  that  the  imaginary  part  of 
the  pole  p”  increases  when  the  electrical  thickness  or  the  loss  tangent  of  the  substrate  is 


29 


increased.  However,  for  a  realistic  case  where  d/Xo<0.1  and  tan5=10'4,  p"  is  small  and 
the  pole  remains  essentially  real.  In  this  case,  the  Green's  function  is  still  sharply  peaked 
near  the  singularity.  In  Figs.  4.2-4.4,  the  component  GXx  of  the  Green’s  function  is 
plotted  for  values  of  p  near  the  singularity.  Even  for  the  thickest  substrate  in  Fig.  4.4,  the 
Green’s  function  is  not  smooth  enough  to  numerically  integrate  without  an  adaptive 
routine.  Therefore,  adding  loss  to  the  dielectric  is  not  a  viable  technique. 

Instead  of  integrating  along  the  real  p  axis,  an  alternate  integration  path  might  be 
taken  around  the  pole  by  applying  the  Cauchy-Goursat  theorem  [10].  This  theorem  states 
that  if  a  function  f(z)  is  analytic  at  all  points  interior  to  and  on  a  simple  closed  contour  C, 
then  the  contour  integration  of  f(z)  on  C  is  zero.  Expressed  in  another  way,  a  line  integral 
of  f(z)  is  independent  of  path.  In  Fig.  4.5,  the  integration  path  b,c,d  in  the  complex  p 
plane  is  then  equivalent  to  the  original  path  a,  provided  the  integrand  is  analytic  and  the 
dielectric  is  slightly  lossy. 

When  this  procedure  was  implemented  in  the  evaluation  of  the  inner  products,  it 
was  found  that  the  results  obtained  were  incorrect.  For  instance,  all  of  the  inner  products 
in  which  the  testing  function  is  identical  to  the  basis  function  should  be  equal.  However, 
all  of  these  inner  products,  known  as  self-terms,  were  different  when  the  integration  path 
was  deformed. 

The  problem  with  the  contour  deformation  is  better  understood  by  first  examining 
the  moment  method  procedure  in  the  spatial  domain,  and  then  determining  how  this  relates 
to  the  spectral  domain  approach.  In  the  spatial  domain,  the  inner  product  is  defined  without 
the  conjugation.  As  an  example,  the  inner  product  for  the  x  component  is 

oo  oo 

(jx»Ex)  =  J  Jjx(x,z)Ex(x,z)  dxdz 


(4.12) 
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where  E*  is  the  electric  field  due  to  the  current  source  a.  The  inner  product  integral  in 
(4.12)  is  equivalent  to  a  reaction  [1 1].  The  reaction  concept  is  useful,  since  it  is  consistent 
with  the  reciprocity  theorem.  For  instance,  the  reaction  of  field  a  on  source  b  is  equal  to  the 
reaction  of  field  b  on  source  a.  The  resulting  symmetries  in  the  testing  procedure  can  be 
used  to  reduce  the  number  of  inner  products  that  need  to  be  evaluated. 

The  physical  insight  gained  from  the  spatial  domain  approach  is  directly  applicable 
in  the  spectral  domain,  because  these  two  approaches  are  actually  equivalent.  This 
equivalence  is  shown  by  applying  Parseval's  theorem  to  Equation  (4.12)  to  give 

oo  oo 

(J5,E*)  =  --p- J  pr(a,p)E|(a,p)docd(3  (4.14) 

— oo  —  oo 

where 

E^(a,P)  =  Gxx(a,p)j^(a,p)  (4.15) 

Here,  the  spatial  domain  testing  function  J*  was  assumed  to  be  real,  and  the  conjugation  of 
T*  is  a  result  of  Parseval's  theorem.  Except  for  the  constant,  the  testing  process  in  (4. 14)  is 
the  same  as  for  the  inner  product  definition  (3.7)  adopted  in  this  work.  Therefore,  the 
spectral  and  spatial  domain  approaches  are  analogous  with  this  inner  product  definition. 

Unfortunately,  if  the  inner  product  definition  is  taken  exactly  as  given  in  (3.7),  it 
can  be  shown  that  the  integrand  is  not  analytic.  The  problem  arises  from  the  conjugation  of 
the  testing  function  in  the  inner  product  definition.  When  the  integration  variables  are 
considered  complex,  the  conjugated  testing  functions  fail  to  satisfy  the  Cauchy-Riemann 
equations  [12],  a  necessary  condition  for  a  function  to  be  analytic.  Since  the  integrand  is 
not  analytic,  the  contour  deformation  technique  cannot  be  applied. 

Although  the  contour  deformation  technique  might  be  valid  with  a  different  inner 
product  definition,  we  chose  to  integrate  the  singularity  with  an  adaptive  routine.  In  this 
way,  the  relationship  between  the  spectral  domain  approach  and  the  spatial  domain 


approach  is  definitely  preserved.  The  numerical  integrator  used  in  the  region  near  the 
singularity  was  an  adaptive  quadrature  routine  designed  specifically  for  singularities. 
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4.3  Convergence  Rate 

When  the  basis  functions  described  in  Chapter  3  are  used,  the  convergence  rate  of 
the  inner  product  integrals  is  extremely  slow.  The  primary  reason  for  this  is  that  the 

integrations  over  the  p  variable  converge  slowly  near  the  a  and  (3  axes.  Each  of  the  basis 

—  —  *1 
functions  Jx,  and  J2j  in  Equations  (3.16)  and  (3.17)  contain  a  p  term  in  polar  coordinates, 

but  their  actual  rate  of  decay  is  much  less  than  p  for  some  0.  For  instance,  the  triangle 

function  in  Jxi  will  tend  to  a  constant  near  0=0'  (the  (3  axis)  resulting  in  a  rate  of  decay  only 

slightly  greater  than  p‘l.  Another  reason  for  the  slow  rate  of  convergence  is  that  the  basis 

functions  are  typically  small  in  the  spatial  domain.  The  corresponding  spectral  domain 

basis  functions  are  very  broad  and  slowly  decaying. 

In  an  effort  to  increase  the  rate  of  convergence,  higher-order  basis  functions  were 
considered.  The  pulse  function  of  (3.12)  was  replaced  by  a  pulse  convolved  with  a 
triangle,  and  the  triangle  function  of  (3.13)  was  replaced  by  a  triangle  convolved  with  a 
triangle.  Each  of  the  resulting  basis  functions  contributed  a  p'7  term  to  the  integrand,  and 
the  inner  product  computation  time  was  reduced. 

However,  this  method  of  increasing  the  rate  of  convergence  suffers  from  the  same 
problems  associated  with  the  original  basis  functions.  The  integrations  over  the  p  variable 
still  converge  more  slowly  near  the  a  and  [3  axes,  and  inner  products  with  small  basis 
functions  tend  to  converge  slowly.  Since  the  rate  of  convergence  with  higher-order  basis 
functions  is  still  inadequate,  a  different  method  of  accelerating  the  integrations  is  discussed 
in  detail  in  Chapter  5.  This  method  is  based  on  extracting  the  asymptotic  behavior  of  the 
integrand. 
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Figure  4. 1 .  Location  of  the  first  TM  surface  wave  pole  for  various  substrate  thicknesses 
and  loss  tangents  (£r=9.6,  d=25  mils). 
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Figure  4.2.  Magnitude  of  the  Green's  function  near  the  TM  surface  wave  pole  for 
dAo=0.02  (er=9.6,  d=25  mils,  9=45*). 
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Figure  4.3.  Magnitude  of  the  Green's  function  near  the  TM  surface  wave  pole  for 
dAo=0.04  (er=9.6,  d=25  mils,  0=45'). 
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Figure  4.4.  Magnitude  of  the  Green's  function  near  the  TM  surface  wave  pole  for 
dAo=0.08  (er=9.6,  d=25  mils,  9=45'). 
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Figure  4.5.  Original  integration  path  a,  through  the  singularity,  and  the  alternate 
integration  path  b,c,d  in  the  complex  p  plane. 
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CHAPTER  5 

ASYMPTOTIC  ACCELERATION  TECHNIQUE 

The  most  serious  limitation  of  the  spectral  domain  approach  is  the  computation  time 
required  to  evaluate  the  inner  products.  As  discussed  earlier,  the  inner  products  are  doubly 
infinite  integrals  that  are  usually  very  slowly  convergent.  The  convergence  problem  is  not 
unique  to  the  analysis  of  microstrip  discontinuities,  and  work  has  been  done  to  reduce  the 
computation  time  in  the  analysis  of  periodic  problems  [2],  such  as  arrays  of  scatterers  and 
arrays  of  antennas.  Although  the  inner  products  in  these  periodic  problems  are  infinite 
sums  rather  than  integrals,  the  general  idea  is  the  same.  A  suitable  function  is  first 
subtracted  from  the  integrand  so  that  the  integral  will  converge  more  rapidly.  The  function 
is  then  integrated  separately  and  added  back  so  that  the  original  integral  remains  unchanged. 

In  this  chapter,  the  acceleration  technique  is  applied  to  the  inner  product  integrals  in 
order  to  increase  their  rate  of  convergence  for  the  infinite  rho  integrations.  In  Section  5.1, 
the  asymptotic  form  of  the  integrand  is  derived  and  then  subtracted  from  the  original 
integrand.  As  a  consequence,  the  resulting  integral  will  converge  more  rapidly.  In  order  to 
recover  the  original  inner  product,  the  asymptotic  form  is  integrated  separately  and  added 
back.  In  Section  5.2,  the  integral  of  the  asymptotic  form  is  evaluated  in  closed  form  by 
deriving  an  equivalent  integral  in  the  spatial  domain. 

5.1  Asymptotic  Form  in  the  Spectral  Domain 

In  applying  the  asymptotic  acceleration  technique,  it  is  important  to  realize  that  any 
asymptotic  expression  can  be  subtracted  from  the  integrand  and  added  back  on,  without 
changing  the  original  integral.  However,  a  speedup  in  the  overall  computation  is  achieved 
only  if  the  expression  that  is  added  back  can  be  integrated  efficiently.  With  this  in  mind, 
the  asymptotic  form  of  the  integrand  is  derived  by  determining  the  asymptotic  behavior  of 
the  Green's  function  and  then  multiplying  this  by  the  basis  and  testing  functions.  In  this 


way,  the  asymptotic  form  is  simplified,  and  the  integral  of  the  asymptotic  form  may  be 
evaluated  in  closed  form. 
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The  derivation  of  the  asymptotic  form  of  the  Green's  function  is  relatively  straight¬ 
forward.  The  behavior  of  each  term  of  the  Green's  function  in  (2.6)-(2.9)  is  determined 
first,  and  then  the  overall  asymptotic  behavior  of  the  Green's  function  is  deduced.  For 
large  p,  the  terms  s  and  s'  are,  to  first  order, 


2  =  -Wp2“ko  =-jP  (5.1) 

s'  =  ~Wp2  ~ko  *  “jP  (5-2) 

Inserting  (5.1)  and  (5.2)  into  the  terms  Dj,  D3  and  D3  yields 

Di=-jp(l  +  er)  (5.3) 

D2  =  — 2  jp  (5.4) 

D3  =  -2jp  (5.5) 


With  the  substitution,  the  term  cot(s'd)  in  Dj,  D2,  and  D3  becomes  jcoth(pd).  The 
hyperbolic  cotangent  is  then  approximated  as  1  for  large  p.  Finally,  Equations  (5.3)-(5.4) 
are  substituted  into  the  expressions  for  the  Green's  function  in  (2.6)-(2.9),  resulting  in  the 
asymptotic  Green's  function  denoted  as 


GxX(P-e)  = 


_  J 


2-2, 


coe0p 


p  sin  (8)  kg 
1  +  8r  2 


2 


Gxz(P>9)  = 


_  J 


2 


OJ£0p 


p^  sin(0)cos(9) 
l  +  er 


(5.6) 

(5.7) 


G^x(p,e)=G^(p,9) 


(5.8) 
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Gzz(p.e)  = 


j 

coe0p 


^p2  cos2  (9) 
1  +  £r 


(5.9) 


The  asymptotic  form  of  the  integrand  is  derived  by  multiplying  the  asymptotic 
Green's  function  in  (5.6)-(5.9)  by  the  appropriate  basis  and  testing  functions  in  the  inner 
product.  This  asymptotic  form  is  then  subtracted  from  the  original  inner  product  integrals 
and  added  back  to  give 


2tz°° 

^xj’GxxJxi  i 

^  =  J  J  ^xj(Gxx~GXx 

)jxi  pdpd0  +  J  |  JxjG^x  Jxi  pdpde 

(5.10) 

0  0 

0  0 

2tc°° 

2?t=« 

/xj’Gxz^zi) 

=  J  f  JljfGxz-Gfx 

)  J zi  pdpd6+  J  j"  jxjGXz  ^zi  pdpdG 

(5.11) 

0  0 

0  0 

27C°o 

2te°o 

^zj’Gzz-^zi) 

=  1 J  Ili(Ozz-GL) 

jzi  pdpde  +  J  |  Jzj Gzz  Jzi  pdpde 

(5.12) 

0  0 

0  0 

where  the  arguments  of  the  functions  have  been  dropped  for  brevity.  The  inner  product 
<JZJ,GZXJX1>  will  no  longer  be  discussed,  since  it  is  identical  in  form  to  (5.1 1).  Although 
the  limits  of  integration  shown  in  (5.10)-(5.12)  encompass  the  entire  plane,  these  limits  can 
be  reduced  to  only  the  first  quadrant.  The  results  of  Section  4.1  can  be  applied  directly  to 
the  integrals,  since  the  symmetries  in  the  integrands  are  unchanged. 

By  subtracting  the  asymptotic  form,  the  rate  of  convergence  for  the  first  integral  in 
(5.10M5.12)  is  accelerated.  As  the  Green's  function  approaches  its  asymptotic  limit  for 
large  p,  the  entire  integrand  will  go  to  zero.  The  second  integral  in  each  inner  product  is 
the  asymptotic  form  that  has  been  added  back.  In  order  to  substantially  reduce  the  overall 
computation  time,  the  second  integral  must  be  computed  efficiently.  In  the  following 
section,  these  integrals  are  evaluated  in  closed  form  in  the  spatial  domain. 

It  is  important  to  note  that  the  asymptotic  form  derived  in  this  section  is  a  first-order 
expression.  For  instance,  the  radicals  in  (5.1)  and  (5.2)  and  the  term  coth(pd)  were 
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approximated  by  only  the  first  term  in  the  series  expansions  for  these  functions.  Although 
an  asymptotic  form  containing  higher-order  terms  will  increase  the  rate  of  convergence  of 
the  first  integral  in  (5.10)-(5.12),  a  price  is  paid  since  the  second  integral  will  be 
considerably  more  complicated.  If  the  second  integral  must  be  evaluated  numerically,  little 
or  no  speedup  in  the  overall  computation  may  be  realized. 

Even  with  a  first-order  approximation,  the  Green's  function  approaches  its 
asymptotic  form  fairly  rapidly.  As  an  example,  the  component  Gxx  of  the  Green's  function 
is  plotted  along  with  Gxx-Gxx  for  different  theta  cuts  in  Figs.  5. 1-5.3.  It  is  apparent  from 
these  plots  that  the  asymptotic  acceleration  technique  will  increase  the  rate  of  convergence 
for  the  rho  integrations  for  all  values  of  theta.  This  is  a  substantial  improvement  over  the 
technique  discussed  in  Section  4.3  in  which  higher-order  basis  functions  are  used  to 
accelerate  the  rate  of  convergence.  With  that  technique,  the  rate  of  convergence  was  not 
improved  significantly  for  values  of  theta  near  0*  or  90*.  Another  advantage  of  the 
asymptotic  acceleration  technique  is  that  the  acceleration  is  dependent  on  the  Green's 
function  and  not  on  the  basis  functions.  Therefore,  the  inner  product  integrations  with 
small  spatial  domain  basis  functions  will  also  converge  rapidly,  even  though  the 
corresponding  spectral  domain  basis  functions  are  very  broad. 

5.2  Analytic  Integration  of  the  Asymptotic  Form 

The  most  difficult  step  in  applying  the  asymptotic  acceleration  technique  is 
evaluating  the  integral  of  the  asymptotic  form.  This  integral  is  the  second  integral  in 
Equations  (5.10)-(5.12).  The  procedure  used  to  evaluate  the  integrals  is  summarized  as 
follows.  The  terms  in  the  integrand  are  first  grouped  into  two  functions,  and  then 
Parseval's  theorem  is  applied  to  derive  an  equivalent  integral  in  the  spatial  domain.  The 
advantage  of  evaluating  the  integral  in  the  spatial  domain  is  that  the  limits  of  integration  are 
then  finite.  At  this  point,  the  integrals  might  be  computed  numerically.  However,  these 
spatial  domain  integrals  are  evaluated  analytically  in  closed  form.  In  this  way,  the 


computation  time  is  reduced  to  a  minimum  and  no  additional  error  is  introduced  in  adding 
back  the  asymptotic  form. 

In  order  to  develop  the  notation  for  the  following  discussion,  Parseval's  theorem  is 
first  stated:  If  U(a,P)  and  V(<x,P)  are  Fourier  transforms  of  u(x,z)  and  v(x,z), 
respectively,  then 


OO  OO  CO  oo 

I J  u*(x,z)v(x,z)  dxdz  ~~^~2  J 1  U*(cc,P)V(a,p)  dad(3  (5.13) 


— oo  — oo 


— oo  ~oo 


In  the  analysis,  the  integrals  in  (5.10)-(5.12)  will  be  considered  in  Cartesian  coordinates  so 
that  Parseval's  theorem  and  other  Fourier  transform  properties  may  be  applied. 

The  terms  in  each  integrand  are  now  grouped  to  form  two  functions  corresponding 
to  the  functions  U*(a,P)  and  V(a,P)  in  (5.13).  When  grouping  the  terms,  particular 
attention  must  be  paid  to  insure  that  each  of  the  resulting  functions  has  an  inverse  transform 
that  exists.  With  this  in  mind,  the  asymptotic  integrals  in  (5.10M5.12)  are  rewritten  as 


jxj.Gxxjxi)=  J  J  T^=  (a2 A2Wx(a)n?x(P)e-^'+Pz'>)dodp 


2mo 


J  J  -r==j  (A2Wx(a)n2x(P)e-J(ax/+Pz/))dadp 


(5.14) 


T  Ga  J 
J  xj  > VJ  xz  J  zi 


OO  oo 

\-  J  [  f 

f  \ 

1 

'  aieo(l  +  er)_J00-oo 

(apAWx(a)fltz(a) 
x  Aw(P)f[t  (P)e-j(ax,+Pz,))dotdP 


(5.15) 
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.jeL 

2ooen 


oo  oo 

J  J 

— oo  — oo 

f  \ 

1 

Jcc2+p2y 

oo  oo 

s  s 

— oo  — oo 

(  \ 

l 

( w-> - j  1  ft,® 

(a12  (p)  n?z(a)  e--i(ax'+Pz,)  )dadp 


(5.16) 


where 


x'  =  (Xi  -Xj) 
z'  =  (Zi-Zj) 


and  where  the  functions  corresponding  to  U*(a,P)  and  V(a,P)  in  (5.13)  are  separated  by 
parentheses.  The  terms  a2,P2,  and  aP  from  the  asymptotic  Green's  function  have  been 
grouped  with  the  basis  and  testing  functions  to  form  the  function  V(a,P).  The  terms  were 
grouped  in  this  manner,  because  the  inverse  transform  of  the  entire  asymptotic  Green's 
function  does  not  exist.  The  remaining  function  1 W  a2+p2  corresponds  to  the  function 
C*(a,P).  Since  this  function  is  real,  the  conjugation  is  dropped  to  simplify  the  notation. 

The  functions  u(x,z)  and  v(x,z),  the  inverse  transforms  of  the  functions  U(a,p) 
and  V(a,P),  are  needed  to  continue  the  derivation  of  the  equivalent  spatial  domain  integral. 
The  inverse  transform  of  U  is  considered  first  and  is  given  by 


U<x.z)  •-K]]  T  1  ,  e^^dadp 

4lt  J_-_i/a2+P2 


(5.17) 


The  integral  in  (5.17)  is  then  converted  into  polar  coordinates  in  both  the  spectral  and 
spatial  domain  variables  [13].  Making  the  substitutions  a=psin9,  P=pcos9,  x=rsinp,  and 
z=rcos<t)  leads  to 


u(r,<J>)  = 


27t°° 

II 


0  0 


ejprcos(e-4»dpd0 


(5.18) 
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where  a  trigonometric  identity  was  used  to  simplify  the  argument  of  the  exponent.  The 
integral  over  the  variable  0  in  (5.18)  is  recognized  as  the  Bessel  function  J0(pr).  Making 

oo 

this  substitution  and  using  the  fact  that  Jj0(t)dt=l  yield  the  inverse  transform 

o 

u(x,z)  = - r-!,— ---  (5.19) 

27TVX2  +Z2 

Determining  the  inverse  transforms  of  the  functions  corresponding  to  V(a,P)  is 
considerably  more  tedious.  In  the  spectral  domain,  these  functions  are  denoted  as 


Vi(a,P)  =  a2  A2Wx(a)n2x(P)e-j(ax'+l5z,) 

(5.20) 

V2(a,P)  =  A2Wx(a)n2x(P)e-j(ax'+f3z') 

(5.21) 

V3(a,P)  =  aPAWx(a)ntz(a)AWz(p)ntx(p)e-j(ax'+Pz,> 

(5.22) 

V4(a,P)  =  p2  Aiz(P)fl2  (a)e-J(ax'+pz,) 

(5.23) 

V5(a,p)  =  A2Wz(P)fl2  (a)e-j(ax'+Pz'> 

(5.24) 

After  applying  several  Fourier  transform  properties,  the  inverse  transforms  of  the  functions 
in  Equations  (5  20)-(5.24)  can  be  represented  as 


^2  p  - 

vj(x,z)  =  “^j[Awx (x  ~  x')* Aw, (x  -  x')  ntx(z-z')*ntx(z-z/)j  (5.25) 

v2(x,z)  =[aWx(x-x,)*Awx(x-x')  ntx(z-z>ntx(z-z')]  (5.26) 

a2  r  1 

v3(x,z)  =  -^lAwx(x  - x')*ntz(x  -  x')  AWz(z-  z>ntx(z  -  z')J  (5.27) 

a2  r  T 

v4(x,z)  =  ~^T[ntz(x -  x')*ntz(x  - x')  AWz(z  -  z> AWz(z  -  z')J  (5.28) 

v5(x,z)  =  [ntz(x-x')*ntz(x-x')  AWz(z-z')*AWz(z-z')]  (5.29) 
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In  (5.25M5.26),  the  partial  derivatives  in  the  spatial  domain  correspond  to  the 
multiplication  by  a  or  (5  in  the  spectral  domain.  The  convolutions  of  pulse  and  triangle 
functions  result  from  the  multiplication  of  their  transforms  in  the  spectral  domain.  Each  of 
these  convolutions  is  performed  over  only  one  spatial  variable,  because  the  functions  are 
separable.  Finally,  the  arguments  of  these  functions  are  shifted  in  the  spatial  domain  due  to 
the  exponential  term  in  the  spectral  domain. 

At  this  point,  deriving  the  final  expressions  for  the  inverse  transforms  is  primarily 
an  exercise  in  algebra.  The  various  convolutions  between  the  triangle  and  pulse  functions 
that  are  required  are  given  in  Table  5.1.  Inserting  these  convolutions  into  equations  (5.25)- 
(5.29)  and  performing  the  indicated  differentiations  will  lead  to  the  final  expressions  for  the 
inverse  transforms  of  v,(x,z).  Since  these  expressions  are  fairly  lengthy,  they  are  not 
presented  here  explicitly. 

Once  the  functions  corresponding  to  u(x,z)  and  v(x,z)  are  determined,  the 
equivalent  asymptotic  integrals  in  the  spatial  domain  are  known  from  (5.13).  These  spatial 
domain  integrals  have  finite  limits  of  the  integration  that  are  determined  by  the  convolutions 
of  the  pulse  and  triangle  functions.  Since  these  convolutions  result  in  piecewise  continuous 
functions,  the  integrals  are  actually  broken  up  into  a  sum  of  integrals  with  limits  over 
smaller  rectangular  regions. 

The  types  of  integrals  that  need  to  be  evaluated  are  polynomials  in  x  and  z  divided 
by  the  function  1/V  x2+z2.  The  polynomials  are  a  result  of  the  convolutions  in  the 
functions  v;,  and  the  radical  comes  from  the  function  u  in  Equation  (5.19).  Most  of  these 
integrals  can  be  found  in  standard  integral  tables,  but  several  of  these  integrals  required 
considerable  effort  to  evaluate.  All  of  the  general  integral  types  needed  are  listed  in 
Table  5.2  in  the  form  of  definite  integrals.  In  this  table,  Ii(x1,x2,z1,z2)  denotes  the  ith 
integral  type,  where  (x!,x2)  and  (zj,z2)  are  the  limits  of  the  x  and  z  integrations. 

When  Equations  (5.19)  and  (5.25)-(5.29)  are  combined  together  with  (5.13)  and 
the  integrals  in  Tables  5.1  and  5.2,  a  closed-form  expression  for  the  integral  of  the 
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asymptotic  form  is  obtained.  The  asymptotic  integral  for  the  inner  product  <Jxj,GxxJxi>  is 
expressed  in  compact  form  as 


.  J 

'  coe0(l  +  er) 


l  jkg 

2C0£o 


(5.30) 


where  the  integrals  I  and  I  are  given  in  Table  5.3.  The  expression  for  the  asymptotic 
integral  <Jzj,Gyz,>  is  obtained  from  (5.30)  by  replacing  x  by  z  and  z  by  x.  Finally,  the 
asymptotic  integral  for  the  inner  product  <JXJ,Gx2Jzr>  is 


=  4tr 


cue, 


i(1+er) 


(5.31) 


where  the  integral  I  is  given  in  Table  5.4.  For  this  integral,  there  are  six  different  cases 
depending  on  the  relative  widths  and  thicknesses  of  the  basis  and  testing  functions. 

Since  the  derivation  of  these  expressions  was  rather  involved,  the  final  expressions 
for  the  asymptotic  integrals  were  verified  by  direct  numerical  evaluation  of  the  spectral 
domain  integrals  in  Equations  (5.14)  and  (5.15).  As  a  test  case,  the  substrate  was  chosen 
to  be  1.57  mm  thick  with  a  relative  dielectric  constant  of  4.  The  operating  frequency  is 
1GHz,  and  the  relative  position  (x',z')  of  the  basis  and  testing  functions  is  (0. 3,0.4).  For 
this  test  case,  the  results  obtained  from  the  analytic  expressions  and  the  numerical 
integrations  are  shown  in  Table  5.5.  All  of  the  analytic  and  numerical  results  are  seen  to 
agree  within  2  percent.  The  difference  is  easily  accounted  for  by  the  error  in  the  numerical 
integration,  because  the  integrand  of  the  asymptotic  form  is  poorly  behaved.  Therefore,  the 
analytic  expressions  for  the  asymptotic  integrals  appear  to  be  correct. 

In  Chapter  6,  the  performance  of  the  asymptotic  acceleration  technique  is  evaluated 
by  comparing  the  results  of  a  routine  that  uses  the  technique  and  one  that  does  not. 
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Table  5. 1 .  Various  convolutions  of  the  pulse  function  and  the  triangle  function. 


Convolution  of  a  Pulse  with  a  Pulse: 

T  <  t  <  0 
0  <  t  <  T 


nT(t)*nT(t)  - 


fT  +  t, 
iT-t, 


Convolution  of  a  Triangle  with  a  Triangle: 


Aw(t)*Aw(t) 


(2w  +  t)3 
6w" 

_ £ _ r_+2w 

2w2  w  3 
t3  t2  |  2w 
2w2  w  3 
(2w  - 1)3 
.  6w2  ’ 


-2w  <  t  <  -w 


-w  <  t  <  0 

0  <  t  <  w 

w  <  t  <  2w 


Convolution  of  a  Triangle  with  a  Pulse: 
Case  1:  T<w 


Aw(t)*nT(t) 


(T  +  2w  +  2t)‘ 


8w 


T 


T 

+  — t, 


w 


4w  w  ’ 


(T  +  2w-2t)2 


-(w  +  T/2)  <t<-(w-T/2) 
-(w-T/2)  <  t  <  -T/2 

-T/2<t<T/2 
T/2  <  t  <  (w-T/2) 


8w 


(w-T/2)  <  t  <  (w  +  T/2) 


Table  5.1  (Corn.) 


Case  2:  w<T<2w 


Aw(t)*nx(t)  -  •< 


Case  3:  T>2w 


Aw(t)*nT(t)  -  < 


(T  +  2w  +  2t)“ 

8w 

T  +  w  (T/2  +  t)2 
- + 1  -  ■ 


2w 

^  7  t" 

T  i  ’ 

4w  w 

T  +  w  (T/2-t)2 

2  2w 

(T  +  2w  -2t)2 


8w 


-(w  +  T/2)<t<-T/2 

-T  /  2  <  t  <  -(w  -  T  /  2) 
-(w-T/2)  <t<  (w-T/2) 
(w-T/2)<t<T/2 

T/2<t<(w  +  T/2) 


(T  +  2w  +  2tr 


8w 

T  +  w  (T/2  +  t)2 
- +  t  -  - 


2w 


w, 

T  +  w _ (T  /  2  —  t)2 

2  1  2w 

(T  +  2w  -2t)2 


-(w  +  T/2)<t<-T/2 

-T  /  2  <  t  <  -(T  /  2  -  w) 
-(T/2-w)£t<(T/2-w) 

(T/2-w)<t<T/2 

T/2<t<(w  +  T/2) 
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Table  5.2  Definite  integrals  needed  in  the  integration  of  the  asymptotic  form. 


z2  fx2  1 


r Zi  fX2 

Il(xl.X2,2l,Z2)  =  j2i  jxi-^ 


X2+Z2 


^dxdz 


=  z 


ln^x  +  Vx2  +  z2  j  - 1  +  xln(^z  + Vx2  +  z2 


z2 


J'z2  fx2  X 

z  Jx  Tr=Tdxdz 


zi  Jxi 


zVx2  +  z2  +  x2  Inf  z  +  Vx2  +  z2 


)] 


z2 


i3(x1,x2,zbz2)=j;2j;2^i_.dxd2 


I X  +  z 


3  3 


=  — Vx2  +  z2  -  —  +  —  InfzWx2  +  z2l~—  lnfx  +  Vx2  +z2 
6  9  3  V  ;  6  V 


Z2 


fz2  fx2  x* 

I4(x1,X2,z1,z2)  =  JZ]  Jxi  -j 


2  2 
X  +Z 


-dxdz 


=  (x2  +z2)2  +  — Vx2  +z2  +— lnfz  +  Vx2  +  z2 


z2 


rz2  rx2  xz 

I5(xi,x2>zi,z2)  =  I  I  7  dxdz 

Jzi  Jxi  Vx2+z2 


3 

x2 

(x2+z2f 

X1 

z2 


Z1 
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Table  5.2  (Cont.) 


I6(xifX2,zi,z2)  =  J  2 J 


z2  fx2  X  Z 


Z1  ,/x2  +  z2 


3  „2 


=  -j(x2  +  z2)2  —  'Vx2  +  z2  +  +z2] 


I7(Xi,X2,Z,,Z2)  =  ffff-T-jr  2  dxdZ 


Vx2  +  z2 


=  1(3x2-2z2)(x2  +  z2)^  ^ 
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Table  5.3  Expressions  for  the  integrals  I  and  \  in  Equation  (5.30).  The  definite  integrals 
Ii  appearing  in  this  table  are  defined  in  Table  5.2. 


I  =  — -(a  +  b) 

2rt 

I  =  — (a  +  b  +  c  +  d) 

2  1 

where 

4 

a  =  Xai[I5(xi’xi+l’zl>z2)_I5(xi’xi+l-z2>z3)  +  'cI2(xi-xi+l-zl-z2)+Kl2(xi>xi+l.z2'z3)] 
i=l 
4 

b  =  Xbi[I2(zl>z2’xi«xi+l)-I2(z2’z3-xi-xi+l)  +  Tll(zl-z2>xi’xi-rl)  +  KIl(z2»z3>xi*xi-l)] 
i=l 

4 

a  =  2iai[l7(Xj,xi+1,z1,Z2)-l7(xi,xi+1,Z2,Z3)  +  tl4(xi,xi+1,z1,z2)  +  Kl4(xi,xi+1,Z2,Z3)] 
i=l 
4 

b  =  Zbi[I6(xi’xi+l-zl>z2)“I6(xi>xi+l*z2*z3)  +  'CI3(xi’xi+l*zl*z2)  +  Kl3(xi>xi+l»z2-z3)] 
i=l 
4 

c  =  Xci[l5(Xj,Xi+1,Z1,Z2)  -  l5(xi,Xj+1,Z2,Z3)  +  Xl2(Xj,Xi+1,Z1,Z2)  +  Kl2(Xj,Xi+1,Z2,Z3)] 
1=1 
4 

d  =  2,di[l2(z1,Z2,xi)xi+1)-I2(Z2,Z3>xi,xi+1)  +  xI1(z1,z2,Xi,xi+1)  +  Kl1(z2,Z3,xi,xi^1)] 
i=i 


T  =  tx-Z' 

K  =  tx+z' 

xj  =  x'  -  2wx 

N 

It 

N 

% 

1 

r» 

x2  =  X'"wx 

z2  =  z' 

X3  =  x' 

Z3  =  z'+  t 

x4  =  x'  +  wx 

X5  =  x'  +  2wx 

and  where 


Table  5.3  (Corn.) 
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al  -  2 

bi  =• 

wx 

-3 

b2  = 

a2  =  — 

wx 

a3  =  — a2 

b3  = 

a4  =  ~al 

b4  = 

2 _ x^ 

2 


w 


w. 


-2  3x' 


w 


X  wx 

-2  3x' 

2 


w 


W, 


WY 


W, 


*1  = 


a->  = 


6wx 

-1 

2w2 


a3  =  -a2 

a4  =  -ai 


bi  = 


£  -1  3x' 

b2  =  —  + 


bi  = 


wx  2wx 
-1  3x/ 

2wx 


w, 


~  1  X 

b4  — - 1- 


w,  2w 


-  „  2x'  x,z 

C!  =2 - + - J 

wr  2wf 


2x' 

3x'2 

c2  ~  — 
wx 

,  2x' 

3x'2 

c3  - 

wx 

2w2 

64  =  -2  - 

2x/ 

x'2 

wx 

2wx 

2  ,  4wx  x'2  x'3 

dj  =  — 2x  H - —  4 - - 

3  w„ 


3  w, 

a, =2^ 


'X 

,/3 


x 

w 


X 

/2 


/3 


2w: 


6w: 


2w2 


*  „  ,  4wx  x'2  x'3 

d4  =  2x  + - -  H - !“ - - 

3  wx  6w" 
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Table  5.4  Expression  for  the  integral  I  in  Equation  (5.31).  The  definite  integrals  I; 
appearing  in  this  table  are  defined  in  Table  5.2. 


I  —  ^  2L2,{aiz[ajx  ^5(xj’xj+l,zi,zi+l)  ■*"  ^jx  ^2(zi’zi+l>xj»xj+l)] 


5  5 


27t 


j-li=l 


+  biz[ajx  ^2 (x j« xj-4-l » zi » zi-t-l )  +  bjx  ^l(xj>xj+l’zi>zi+l)]| 


where  the  variables  Xj,z,,ajx,bjx,aiz,  and  bjz  are  listed  below  for  each  of  the  six  different 
cases. 


Cass.1:  tz<wx  and  tx<wz 


X1  =  x'  ~(wx  -t-  tz/2) 
x2  =  x'-(wx  -tz/2) 
x3  =  x'~tz/2 
x4  =  x'  +  tz/2 
x5  =  x'  +  (wx  -tz/2) 
x6  =  x'  +  (wx  +  tz/2) 


alx  = - 


w. 


bix=Tk- 
2  w. 


•  +  1 


w, 


a2x  =0 


b2x=-k 

WT 


a3x  = 


w. 


h  -2X' 
b3x  -  — 

wT 


a4x  =0 


t>4x  =  — 
W, 


a5x  = - 

w„ 


.  -tz  X  . 

b5x  =- —  + - 1 

2wx  wx 


Z1  =  Z'-(wz+tx/2) 
z2  =z'-(w2-tx/2) 
z3  =  z'-tx/2 
z4  =z'  +  *  J2 

z5  =  z'  +  (wz-tx/2) 
z6  =  z'  +  (wz  +  tx/2) 


I 

b  -  tx 

alz  - 

wz 

2wz 

a2z  =0 

b2z  =  is- 

wz 

-2 

a3z  = - 

wz 

h  -2z' 
^3z  ~ 

wz 

o 

II 

NJ 

b.‘2  =  — S- 

wz 

1 

b<  -  _tx 

a5z  “ 

wz 

d5z  ~ 

2wz 
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Table  5.4  (Cont.) 

Case  2:  wx<tz<2wx  and  wz<tx<2wz 

*1  =  x'-(wx+tz/2) 

x2  =  x'  ~  tz/2 

x3  =  x'~(wx-tz/2) 

x4  =  x'  +  (wx  -tz/2) 

x5  =  x'  +  tz/2 

x6  =  x'  +  (wx  +  tz/2) 


Z1  =  z'-(wz  +tx/2) 

z2  =  z'-tx/2 

z3  =  z'-(wz-tx/2) 
z4  =z'  +  (wz-tx/2) 

z5  =  z'  +  tx/2 
z6  =  z'  +  (wz  +  tx/2) 


1 

alx  - 

wx 


a2x  - 


~1 

wx 

-7 


a3x  - - 

wx 

~1 

a4x  — 

W„ 


a5x 


1 


w. 


Hx 


b3x  = 
b4x  = 
b5x  = 


x' 

+  1 

alz 

1 

hi 

.  tx 

z' 

2wx 

wx 

wz 

Dlz 

2w: 

wz 

~CZ  | 

x' 

+  1 

a2z 

-1 

b2z 

_  _tx 

z' 

2wx 

wx 

wz 

2wz 

I 

wz 

2x' 

wx 

a3z 

= 

-2 

wz 

b3x 

_  2z' 
wz 

Jz_  + 

x' 

-1 

a4z 

-1 

h  * 

tx 

4_ 

z 

2wx 

wx 

wz 

°4z 

2wz 

r 

wz 

~lz 

x' 

-1 

a5z 

1 

b5z 

..  -tx 

z' 

2wx 

wx 

wz 

2wz 

wz 

+  1 


+  1 


-1 


-1 


Case  3:  2wx-tz  2wz-cx 


x2  =  x'  ~  tz/2 


x5  =  x'  +  tz/2 


) 

zi  =  z'-(wz+tx/2) 

N 

to 

II 

N 

1 

rt 

X 

to 

) 

z3  =  z'  +  (wz  +  tx/2) 

) 

z4  =  z'  +  (tx/2-wz) 

z5  =  z'  +  t  x/'2 

•) 

z6  =  z'  +  (wz  +  tx/2) 
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Case  4:  a)  wx<tz<2wx  and  tx<wz 

Take  ajx,bjx,Xj  from  Case  2,  and  take  aiz,biz,Zi  from  Case  1 . 
b)  wz<tx<2wz  and  tz<wx 

Take  aJX,bJX,Xj  from  Case  1,  and  take  a;z,biz,z;  from  Case  2. 


Q»$C.5:  a)  2wx<tz  and  tx<wz 

Take  ajx,Ojx,Xj  from  Case  3,  and  take  a^.b^.Zj  from  Case  1. 
b)  2wz<tx  and  tz<wx 

Take  ajx,bjx,Xj  from  Case  1,  and  take  ajz,bjz,Zj  from  Case  3. 


Case  6:  a)  2wx<tz  and  wz<tx<2wz 

Take  ajx,bjx,Xj  from  Case  3,  and  take  ajz,biz,Zj  from  Case  2. 
b)  2wz<tx  and  wx<tz<2wx 

Take  ajx,bjx,Xj  from  Case  2,  and  take  ajz,bjz,Zi  from  Case  3. 


Table  5.5.  Numerical  verification  of  the  analytic  expressions  for  the  integral  of  the 
asymptotic  form.  The  numerical  result  is  calculated  by  numerically  integrating 
the  spectral  domain  integrals  in  (5.14)  and  (5.15).  The  analytic  result  is 
obtained  from  the  expressions  (5.30)  and  (5.31).  The  relative  position  (x',z') 
of  the  basis  and  testing  functions  is  (0.3, 0.4).  (f=l  GHz,  d=1.57  mm,  and 
£,=4.0) 


Inner  Product 

wx 

h 

W7 

tz 

Analytic 

Numerical 

<  J  xj»GxxJ  xi>. 

3xl0-2 

4x10*2 

— 

— 

-7.150x10-2 

-7.148x10-2 

<Jxj,GX7Jzj>. 

Case  1 

X 

I — * 

o 

4x10-2 

5x10-2 

2x10*2 

-3.145X10*4 

-3.202xl0-4 

Case  2 

" 

8x10-2 

- 

4x10-2 

-1.265x10-2 

-1.250x10-3 

Case  3 

II 

11x10-2 

«i 

6x10-2 

-2.625x10-3 

-2.638x10-3 

Case  4 

If 

4x10*2 

if 

4x10-2 

-6.288xl0-4 

-6.268x1  O*4 

Case  5 

-• 

4x10-2 

it 

6x10-2 

-9.425x1  O'4 

-9.434x1 0-4 

Case  6 

-• 

8x10-2 

ii 

6x10*2 

-1.896x10-3 

-1.890x10-3 

56 


CHAPTER  6 
PERFORMANCE 

In  this  chapter,  the  asymptotic  acceleration  technique  discussed  in  Chapter  5  is 
evaluated  by  comparing  the  results  of  two  different  algorithms.  One  algorithm  computes 
the  inner  products  by  direct  numerical  integration,  and  the  other  employs  the  asymptotic 
acceleration  technique  in  the  numerical  integration.  The  direct  integration  algorithm  is  used 
as  the  standard  for  determining  the  accuracy  of  the  acceleration  technique.  It  also  provides 
a  basis  for  a  comparison  of  the  execution  speeds. 

Before  comparing  the  results,  it  is  necessary  to  understand  the  numerical  integration 
routines  used  to  compute  the  inner  product  integrals.  Both  algorithms  perform  the  two- 
dimensional  integrations  in  polar  coordinates,  and  by  taking  advantage  of  the  symmetry  in 
the  integrands,  the  limits  of  integration  have  been  reduced  to  the  first  quadrant.  The 
integration  is  implemented  numerically  by  stepping  in  the  theta  variable  and  integrating  in 
rho  at  each  step  until  convergence.  The  rho  integrations  are  then  summed  up  over  the  theta 
steps  with  a  basic  Gaussian  quadrature  integration  routine.  For  the  rho  integrations,  a 
general  adaptive  quadrature  routine  is  used  everywhere,  except  in  the  region  of  the 
singularity.  In  this  region,  an  adaptive  routine  designed  specifically  for  singular  functions 
is  used. 

Defining  a  good  convergence  criterion  for  numerical  integrations  is  a  difficult  task 
in  general.  In  this  work,  the  criterion  chosen  for  truncating  the  infinite  rho  integrations  is 
based  on  a  relative  convergence  test.  After  each  step  in  the  rho  integration,  the  most  recent 
approximation  to  the  integral  is  compared  with  the  approximation  from  the  previous  step. 
If  the  relative  difference  is  less  than  a  user- specified  tolerance,  the  integration  is  terminated. 
The  actual  procedure  implemented  requires  the  integration  to  satisfy  this  convergence 
criterion  two  consecutive  times  before  terminating. 
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Although  the  asymptotic  acceleration  technique  was  designed  to  increase  the  rate  of 
convergence  of  the  rho  integrations,  the  technique  actually  accomplishes  more  than  this.  In 
some  cases,  the  inner  product  evaluation  is  not  possible  without  the  acceleration  technique. 
The  reason  for  this  is  that  the  integrand  becomes  highly  oscillatory  with  respect  to  theta  for 
large  rho.  Since  the  rho  integrations  are  truncated  sooner  with  the  acceleration  technique, 
this  highly  oscillatory  region  of  the  integrand  may  be  avoided.  It  is  important  to  note  that 
the  poor  asymptotic  behavior  has  been  captured  analytically  in  the  asymptotic  acceleration 
technique.  Therefore,  the  results  obtained  with  this  technique  are  possibly  more  reliable 
than  those  obtained  by  direct  numerical  integration. 

As  an  example  of  the  integrand’s  poor  behavior,  the  integrand  of  the  inner  product 
<JXJ,GxxTxj>  is  plotted  with  respect  to  theta  in  Figures  6.1  and  6.2  for  rho  at  1000  and 
3000.  These  values  of  rho  are  reasonable,  since  most  of  the  rho  integrations  cannot  be 
truncated  before  3000  for  these  inner  products.  From  the  plots  it  is  apparent  that  the 
integrand  becomes  more  oscillatory  as  rho  increases.  The  behavior  of  the  integrand  will 
also  become  worse  as  the  separation  between  the  basis  and  testing  functions  increases.  In 
Figure  6.1,  the  separation  is  approximately  a  quarter  of  a  wavelength  in  free  space.  The 
behavior  of  the  integrand  in  Figure  6.2  is  substantially  worse,  where  the  separation  is  about 
a  half  wavelength.  This  poor  behavior  arises  from  the  complex  exponential  term  associated 
with  the  basis  and  testing  functions.  As  a  result  of  this  term,  the  direct  integration  routine 
often  converges  slowly  in  the  theta  integration.  When  the  separation  is  too  large,  the  direct 
numerical  integration  completely  fails  to  converge  before  numerical  error  dominates. 

Since  the  direct  integration  of  the  inner  products  fails  for  large  separations,  the 
comparison  of  the  algorithms  is  limited  to  inner  products  with  relatively  small  separations. 
When  a  comparison  is  possible,  the  agreement  between  the  direct  and  the  accelerated 
evaluations  of  the  inner  products  is  excellent.  As  an  example,  consider  the  case  in  which 
the  operating  frequency  is  1  GHz  and  the  dielectric  substrate  is  1.57  mm  thick  with  a 
relative  dielectric  constant  of  4.  The  width  and  thickness  of  the  basis  functions  were 
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chosen  to  be  approximately  a  twentieth  of  a  free-space  wavelength.  Several  inner  products 
were  computed  for  this  case,  and  the  results  are  given  in  Table  6.1.  The  results  from  the 
two  algorithms  are  seen  to  agree  very  well  for  the  various  separations  of  the  source  and 
field  points  considered.  For  all  of  the  cases,  the  complex  relative  difference  is  within  3 
percent.  The  largest  separation  listed  is  approximately  a  quarter  wavelength.  Beyond  this 
distance,  the  direct  integration  algorithm  did  not  converge  with  a  reasonable  number  of 
points  in  the  theta  integration. 

Note  that  the  real  parts  of  the  inner  products  in  Table  6.1  are  identical  for  the  two 
algorithms,  because  the  asymptotic  form  only  contributes  to  the  imaginary  part  of  the  inner 
product.  This  brings  up  the  question  of  the  true  accuracy  of  either  algorithm.  Although  the 
results  of  the  two  algorithms  are  consistent,  it  is  possible  that  both  results  may  be  in  error, 
since  the  algorithms  are  based  on  similar  numerical  integration  routines.  However,  from 
the  results  presented,  it  is  safe  to  conclude  that  the  asymptotic  acceleration  technique  does 
not  introduce  any  new  error  in  the  evaluation  of  the  inner  products. 

Finally,  the  asymptotic  acceleration  technique  is  evaluated  in  terms  of  its  effect  on 
reducing  the  computation  time.  For  the  cases  listed  in  Table  6.1,  the  acceleration  technique 
provided  a  speedup  that  varied  between  two  and  four.  However,  the  time  saved  will 
increase  as  the  limit  of  the  direct  integration  is  pushed  with  larger  separation  distances.  The 
acceleration  technique  will  also  be  more  effective  for  thicker  substrates,  since  the 
approximation  for  the  hyperbolic  cotangent  in  the  asymptotic  form  is  then  better. 

The  effectiveness  of  the  acceleration  technique  can  also  be  examined  by  comparing 
the  maximum  rho  reached  before  the  convergence  criterion  is  satisfied.  In  Figure  6.3,  an 
actual  path  of  integration  in  the  first  quadrant  of  the  (p,9)  plane  is  plotted  for  both 
algorithms.  Each  radial  line  represents  a  rho  integration  for  a  particular  theta  step.  For  the 
direct  integration,  there  are  a  few  rho  integrations  that  do  not  converge  until  after  5000.  In 
comparison,  all  of  the  rho  integrations  for  the  accelerated  algorithm  are  seen  to  converge 
before  1800.  For  this  example,  it  appears  as  though  the  acceleration  technique  should 


provide  a  speedup  of  more  than  the  observed  factor  of  two.  The  actual  reduction  in 
computation  time  is  more  modest,  because  a  significant  portion  of  the  computation  time  in 
each  rho  integration  is  spent  integrating  the  singularity  adaptively. 
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Table  6.1  Comparison  of  the  inner  products  computed  using  a  direct  integration 
algorithm  and  an  accelerated  algorithm.  The  widths  and  thicknesses  of  the 
"rooftop"  basis  functions  are  each  0.05  XQ,  and  the  separation  (x',z')  is  given 
in  terms  of  wavelengths  (d=1.57  mm,  ^=4.0,  f=l  GHz). 
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(b)  Magnitude  of  the  integrand  with  respect  to  9  at  p=3000. 

Figure  6. 1  Behavior  of  the  integrand  for  <Jxj.GxxJxj>  with  separation  variables  x'  and  z 
equal  to  0.2  X0  (d=1.57  mm,  er=4,  f=l  GHz,wx=tx=0.05  X0). 
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(b)  Magnitude  of  the  integrand  with  respect  to  9  at  o=3000. 


Figure  6.2  Behavior  of  the  integrand  for  ^xj^xx^xP*  with  separation  variables  x'  and 
equal  to  0.4  X0  (d=1.57  mm,  6^=4,  f=l  GHz,wx=tx=0.05  il0). 


(b)  Accelerated  algorithm. 


Figure  6.3 


Integration  path  in  the  first  quadrant  of  the  (p,0)  plane  for  the  inner  product 
<Jxj,GxxJxi>  with  separation  variables  x'  and  z  equal  to  0.2  X^  (d=1.57  mm. 
Er=4,  f=l  GHz,wx=tx=0.05  X0). 
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CHAPTER  7 
CONCLUSIONS 

In  this  thesis  the  general  microstrip  discontinuity  problem  was  formulated  in  the 
spectral  domain.  The  moment  method  was  then  used  to  generate  a  set  of  algebraic 
equations  that  could  be  solved  numerically.  When  this  procedure  was  implemented,  it  was 
found  that  there  were  numerical  difficulties  associated  with  the  evaluation  of  the  inner 
products.  One  such  difficulty  was  the  numerical  integration  of  the  singularities 
corresponding  to  the  surface  wave  poles.  The  nature  of  these  singularities  and  possible 
methods  to  integrate  them  were  discussed.  Another  difficulty  encountered  was  the  slow 
convergence  rate  of  the  inner  product  integrals.  The  convergence  rate  of  these  integrals 
was  improved  by  applying  an  asymptotic  acceleration  technique.  In  this  work,  the  integral 
of  the  asymptotic  form  was  derived  analytically  in  closed  form;  therefore,  the  maximum 
speedup  was  achieved  for  the  asymptotic  extraction.  Also,  no  additional  error  was 
introduced  in  the  computation. 

In  order  to  evaluate  the  effectiveness  of  the  asymptotic  acceleration  technique,  two 
different  algorithms  were  developed.  One  algorithm  employed  the  acceleration  technique  in 
the  evaluation  of  the  inner  products,  and  the  other  computed  the  inner  products  by  direct 
numerical  integration.  Both  routines  took  advantage  of  symmetry  in  the  integrands  to 
obtain  an  initial  speedup.  The  results  from  the  two  algorithms  were  in  close  agreement, 
confirming  the  fact  that  the  acceleration  technique  does  not  compromise  the  accuracy.  For 
the  cases  considered,  the  acceleration  technique  reduced  the  overall  computation  time  by  a 
factor  that  varied  between  two  and  four.  It  was  also  found  that  the  acceleration  technique 
allowed  certain  inner  products  to  be  evaluated  numerically  that  appeared  to  be  unattainable 
with  the  direct  integration  routine. 

Although  the  asymptotic  acceleration  technique  was  developed  for  the  microstrip 
discontinuity  problem,  it  should  be  noted  that  the  technique  is  also  useful  in  the  spectral 


domain  analysis  of  other  grounded  dielectric  slab  structures.  In  particular,  patch  antennas, 
printed  dipoles,  and  disk  resonators  can  be  analyzed  using  this  technique. 
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The  work  to  improve  the  efficiency  of  the  inner  product  evaluation  is  far  from 
complete.  Several  potential  improvements  on  the  technique  should  be  investigated  in  the 
future.  For  instance,  the  acceleration  technique  could  easily  be  implemented  in  conjunction 
with  higher-order  basis  functions  to  further  enhance  the  convergence  rate  of  the  inner 
product  integrals.  In  addition,  the  use  of  higher-order  asymptotic  forms  in  the  acceleration 
technique  should  be  considered.  In  the  latter  case,  an  analytic  evaluation  of  the  resulting 
asymptotic  integrals  may  not  be  possible,  but  approximations  may  exist  (i.e.,  series 
expansions)  for  which  the  tradeoffs  between  computation  time  and  accuracy  are  acceptable. 
Finally,  the  computation  time  could  be  reduced  significantly  if  a  more  efficient  method  of 
integrating  the  singularities  were  implemented. 
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